clear all;
clc;

load mnist0v1;

X = double(data(1:500,:)); %size n, both labeled and unlabeled data
Y = double(label(1:500)); %size l < n, labels

l = size(Y,1);
n = size(X,1);

k = 4;

[IDX,D] = knnsearch(X,X,'K',k);
D2 = D.^2;
IDXvec = reshape(IDX(:,2:k),[],1);
Dvec = reshape(D2(:,2:k),[],1);
Ivec = reshape(repmat((1:l)',1,k-1),[],1);

Dvec = Dvec/max(Dvec);

cvx_begin sdp
cvx_solver sedumi
    variable K(n,n) symmetric
    minimize -trace(K)
    subject to
        sum(sum(K))/n == 0
        K >= 0
        DK = diag(K)
        DKvec = reshape(repmat(DK,1,k-1),[],1)
        DKvec + DK(IDXvec) - 2*K(sub2ind(Ivec,IDXvec)) <= Dvec
cvx_end